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ABSTRACT 



Following Dexter & Agol (2009) we present a new public code for the fast cal- 



culation of null geodesies in the Kerr spacetime. Using Weierstrass' and Jacobi's 
elliptic functions, we express all coordinates and affine parameters as analytical 
and numerical functions of a parameter p ) which is an integral value along the 
geodesic. This is a main difference of our code compares with previous similar 
ones. The advantage of this treatment is that the information about the turning 
points do not need to be specified in advance by the user, and many applications 
such as imaging, the calculation of line profiles or the observer-emitter problem, 
etc become root finding problems. All elliptic integrations are computed by Carl- 



son's elliptic integral method as Dexter & Agol (2009) did, which guarantees the 



fast computational speed of our code. The formulae to compute the constants of 



motion given by Cunningham & Bardeen (1973) have been extended, which al- 



low one readily to handle the situations, in which the emitter or the observer has 
arbitrary distance and motion state with respect to the central compact object. 
The validation of the code has been extensively tested by its application to toy 
problems from the literature. The source FORTRAN code is freely available for 
download on the webQ 
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Introduction 



There are wide interests in calculating the radiative transfer near the compact objects, 
such as black hole, neutron star and white dwarf. The radiation will be affected by various 
effects, such as, light bending or focusing, time dilation, Doppler boosting and gravitational 
redshift, under the strong gravitational field of the compact objects. The consideration of 
these effects not only help us to understand the observed results, therefor to study these 
compact objects, but also even to test the correctness of the general relativity under strong 
gravity. A good example is the study of the fluorescent iron line in the X-ray band at 6.4-6.9 



keV, which is seen in many active galactic nuclei especially for Seyfert galaxies (Laor 1991 



galaxies ( 


Laor 


1991 


Miniutti et al.|2004) 



The line appears broadened and skewed as a result of the Doppler effect and gravitational 
redshift, thus it is an important diagnostic to study the geometry and other properties of the 
accretion flow at the vicinity of the central black hole ( Fabian et al.p OOO). Another example 



is the study of SMBH in the center of our galaxy. It has been comprehensively accepted that 
in the center of our galaxy a super-massive black hole with ~ 4x 10 6 M exists (Schodel et al. 



2003; Gebhardt et al. 2000; Hopkins et al. 2008) and its shadow may be observed directly 



in the radio band in the near future. Based on the general relativistic numerical simulations 



of the accretion flow, Noble et al. (2007) present the first dynamically self-consistent models 



of millimeter and sub-millimeter emission from Sgr A*. Yuan et al. (2009) calculated the 



observed images of Sgr A* with a fully general relativistic radiative inefficient accretion flow. 

A natural way to include all of the gravitational effects is to track the ray following its 
null geodesic orbit. Which requires the fast and accurate computations of the trajectory of 



a photon in the Kerr spacetime. Fanton et al. (1997) proposed a fast code to calculate the 



accreting lines and thin disk images. Cadez et al. (1998) translated all integrations into the 



Legendre's standard elliptic integrals and wrote a fast numerical code. Dexter & Agol (2009) 
presented a new fast public code, named geokerr, for the computing of all coordinates of null 
geodesies in the Kerr spacetime by using the Carlson's elliptic integrals semi-analytically for 
the first time. There are two computational methods often used in these codes, they are: (1) 



the elliptic function method, which relies on the integrability of the geodesies (Cunningham 



& Bardeen 


1973; Cunningham 


1975; Rauch & Blandford 


1994 


ISpeith et al. 1995; I 


7 anton 


et al. 


1997| 


Cadez et al. 1998; : 


A et al. 


2005; Wu & Wang 


2007 


Dexter & Agol 


2009 


Yuan 


et al. 


2009D, and (2) the direct geodesic integration method (Fuerst & Wu||2004; 


Schnittman 


& Rezzolla| |2006; Dolence et al. 2009; 


Anderson et al. 2010; 


Vincent et al. 2011; " 


Y)unsi 



et al. 2012). Usually people regard the direct geodesic integration method to be a better 
choice than the elliptic function method, for the direct geodesic integration method can 
deal with any three-dimensional accretion flow (Younsi et al. 2012| ), especially in radiative 
transfer problems which require the calculations of many points along each geodesic, the 
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direct integration method is simpler and faster (Dolence et al. 2009). While the the elliptic 



function method is considered to be just efficient for the calculation of the emissions come 
from an optically thick, geometrical thin and axisymmetric accretion disk system. But we 



think that the elliptic function method based on Dexter & Agol (2009) after some extensions 



can overcome these shortages and not only handle any three-dimensional accretion flows 
readily, but also be more efficient, flexible, and accurate, because it can compute arbitrary 
points on arbitrary sections for any geodesies. The speed of the code based on this approach is 
still very fast for many potential applications. While the direct geodesic integration method 
must integrate the geodesic from the initial position to the interested points, the waste of 
computational time is inevitable. 

We present a new public code for the computation of null geodesies in the Kerr spacetime 
following the work of Dexter & Agol ( 2009[ ). In our code all coordinates and the affine 
parameters are expressed as functions of a parameter p ) which corresponds to I u or 1^ in 



Dexter & Agol (2009). Using parameter p as the independent variable, the computations 



are easier and simpler, thus more convenient, mainly due to the fact that the information 



about the turning points does not need to be prescribed in advance comparing with Dexter 



& Agol (2009). Meanwhile Yuan et al. (2009) have demonstrated that the parameter p can 



replace the affine parameter to be the independent variable in radiative transfer problems. 
We not only take this replacing, but also used it to handle more sophisticated applications. 
We extend the formulae of computing constants of motion from initial conditions to a more 
general form, which can readily handle the cases in which the emitter or the observer has 
arbitrary motion state and distance with respect to the central black hole. We reduce 
the elliptic integrals from the motion equations derived from the Hamilton- Jacobi equation 
( Carter|1968 ) to the Weierstrass' elliptic integrals rather than to the Legendre's ones, because 
the former ones are much easier to handle. Then we calculate these integrals by Carlson's 
method. 

The paper is organized as follows. In section [2] we give the motion equations for null 
geodesies. In section [3] we express all coordinates and affine parameters as functions of a 
parameter p analytically and numerically. In section [| we give the extended formulae for 
the computation of constants of motion. A brief introduction and discussion about the code 
are given in section [5j In section [6] we demonstrate the testing results of our code for toy 
problems in the literature. The conclusions and discussions are finally presented in section 
[7| The general relativity calculations follow the notational conventions of the text given by 



Misner et al. (1973). The natural unit are used through out this paper, in which constants 



G=c=l, and the mass of the central black hole M is also taken to be 1. 
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2. motion equations 



Following |Bardeen et al.| (1972), we write the Kerr line element in the Boyer-Lindquist 
(B-L) coordinates (£,r, 0,0) as 

ds 2 = -e 2 ^t 2 + e 2 ^(# - cjrft) 2 + e 2 ^dr 2 + e 2 ^ 2 d# 2 , (1) 

where 



e 2 " 
A = r 



EA 



3 2V 



sin" 0A 



A ' " E A' 6 ~ A ' 

2 2r + a 2 , E = r 2 + a 2 cos 2 0, A = (r 2 + a 2 ) 2 - Aa 2 sin 2 0, 



2/^2 



2ar 



(2) 



and a is the spin parameter of the black hole. 

The geodesic equations for a freely test particle read 

d 2 x a _ „ „ 

0, 



d(J 2 + ^ UU 



(3) 



where a is the proper time for particles and an affine parameter for photons, u v is the four- 
velocity, is the connection coefficients. Carter (1968) found these equations are integrable 
under Kerr spacetime and got the differential and integral forms of motion equations for 
particles by using the Hamilton- Jacobi equation. For a photon, whose rest mass m is zero, 
the equations of motion have the following forms: 

_ dr 



da 

E— 
da 

s # 

da 

E— 
da 

± 
a 



— (a 



A , T 
. 9 J + a—, 
sin 2 y A' 



-a(a sin 2 6 - A) + (r 2 + a 2 ) 



T 



dr 



R 



r r 2 



Vr 



t = a + 2 



± 



dr + a 



A\71 



d0 



cos 2 



d0, 



(ir, 



/r rp f>$ 
i=dr + / 
Av^ J 



r " T , ^ Acsc 2 0-a 7/1 
- 7 dr+ I 7= d0, 



(4) 
(5) 
(6) 
(7) 
(8) 

(9) 
(10) 

an 
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where 

R = r 4 - (q + A 2 - a 2 )r 2 + 2[q + (A - a) 2 ]r - a 2 q, (12) 

Qo = q + a 2 cos 2 6 - X 2 cot 2 9, (13) 

T = r 2 + a 2 — Aa, (14) 
q and A are constants of motion denned by 

9= J*' A = ¥' (15) 

where Q is the Carter constant, L z is the angular momentum of the photon about the 
spin axis of the black hole, E is the energy measured by an observer at infinity. The four 
momentum of a photon can be expressed as 



p /1 = J B(-l,±^,± V ^,A), (16) 



which is often used in the discussion of the motion of a photon. 



From the equation (13) we know that if q = and 6 = 7r/2, then Q$ = 0. 6 = 7r/2 means 
the motion of the photon is confined in the equatorial plane forever ( Chandrasekharp 983 ) . 
Thus the motion equations with integral forms now become invalid for Q$ = appeared in 
the denominator. We need the motion equations with differential forms. 



From equation Q we have 

V = / ^T^dr, (17) 



^2 



where the subscript pm means 'plane motion'. Dividing equation ^ by Q and integrating 
both sides, we obtain 



y/R J A ,/r 
Similarly from equation ^ and Q, we obtain 



r A - a , f r T dr 

dr + a -r—i=. (18) 



^ P m = / -^dr + 2 



R J A y/R 

r rT dr 



apm+2 1 ~KTm- (19) 



The spherical motion is an another special case, in which the photon is confined on 
a sphere and the motion of which can be described by equations: R = and dR/dr = 
(Bardeen et al. 1972; Shakura 1987). Thus the motion equations with integral forms also 
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become invalid due to R = appears in the denominator. Similarly from equation ^ we 
have 



dO [ e a 2 cos 2 ,„ 



r d0 r 

J w» + l 



(20) 



where the subscript sm means "spherical motion". Dividing equation ^ by ^ and inte- 
grating both sides, we have 



T 

a— 



a j 

From the equations ^ and ([7]) we have 

4m 0"sm H~ 2r- 



d# A esc 2 0-a in 
+ I — — d0. 



(21) 



T r d0 



A/ 



(22) 



From equation (l8|) , we introduce a new parameter p with following definition to describe 



the motion of a photon along its geodesic (Yuan et al. 2009) 



P 



d0 



(23) 



Because the sign ahead the integral is the same with dr and d6 ) p is always nonnegative and 
increases monotonically as the photon movies along the geodesic. From the above definition, 
we know that r and 9 are functions of p. In the next section, we will give the explicit forms 
of these functions by using Weierstrass' and Jacobi's elliptic functions. 



3. The expressions of all coordinates as functions of p 

3.1. Turning points 

From the equations of motion we know that both R and Qq must be nonnegative. 
This restriction divides the coordinate space into allowed (where R > and Qq > 0) and 
forbidden (where R < or Q$ < 0) regions for the motion of a photon. The boundary 
points of these regions are the so called turning points, their coordinates r tp and 9 tp satisfy 
equations R(r tp ) = and Qo(6 tp ) = 0. For a photon emitted at r ini and 9 ini , its motion will 
be confined between two turning points r tpi and r tp2 for radial coordinate, 9 tpi and 9 tp2 for 
poloidal coordinate. If we assume r tpi < r tp2) and 9 tpi < 9 tp2) then we have r ini G [rt P1 ,r tp2 ] 
and 9i n i G [9t Pl ,9 tp2 \. Because p r = ±V^R/A and pe = ±y / ©^, if p r = (or pe = 0) at 
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the initial position, we have R(r ini ) = (or Qe(9 ini ) = 0), therefore r ini (or 9 ini ) must be a 
turning point and equal to one of r tpi and r tp2 (or 9 tpi and 9 tp2 ). 

The radial motion of a photon can be unbounded, meaning that the photon can go to 
the infinity or fall into the black hole. These cases usually correspond to equation R(r) = 
has no real roots or r tpi is less or equal to the radius of the event horizon. We regard the 
infinity and the event horizon of the black hole as two special turning points in the radial 
motion, a photon will asymptotically approach them but never return from them. Thus r tp2 
can be the infinity and r tpi can be less or equal to (r^ is the radius of the event horizon) . 

For the poloidal motion, there is also two special positions, 9 = and 9 = 7r, i.e., the 
spin axis of the black hole. A photon with A = will go through the spin axis due to zero 
angular momentum, and will change the sign of its angular velocity d9/da instantaneously, 
and its azimuthal coordinate will jump from cj) to cj) ± 7r (Shakura 1987| ), implying that the 



spin axis is not a turning position. From the equation (13) we also know that and tt are 
not the roots of equation Q$ (9) = 0. 



3.2. /x coordinate 



Firstly, we use a new variable ji to replace cos#, and the equation (23) can be rewritten 

as: 

p = ± 



where 

e M = g-( g + A 2 -a 2 V 2 -aV- (25) 

Both R and G M are quartic, but the polynomial of Weierstrass' standard elliptic integral 
is cubic. We need a variable transformation to make R and G M to be cubic. We define the 
following constants for poloidal motion: 

b = -4a 2 /4 pi - 2(q + A 2 - a 2 )^ tpi , (26) 

6 1 = -2^-^ + A 2 -a 2 ), (27) 
4 

b2 = -^a 2 /i tpi , (28) 
h = -a 2 , (29) 
where fi tpi = cos 9 tpi and introduce a new variable t, 

bo 1 bi . . 

4(p-fi tpi ) 4 
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Making transformation from \± to t, the ji part of equation (24) can be reduced to 

rKv) fj+ 
J V 4t - 92t - gs 

where g 2 = \{b\ — 60^2), #3 = ]^(36o&i&2 — 26^ — &0&3). Using the definition of Weierstrass' 
elliptic function jp(z; g 2 , #3) ( Abramowitz fc Stegun|[l965| ), from equation (31), we have t = 
p(p ± 11^; g 2) gz). Solving equation (30) for /x, we can express \i as the function of p: 

Kp) = 777 t tt b ° s — r + ^pi ' ( 32 ) 

where 

n M = |p _1 [t(/i^);^ 2 ,^3]|. (33) 

The sign ahead II M depends on the initial value of po ) which is the 6 component of four 
momentum of a photon, and 



where is the period of <p(z]g 2) g?) and n = 0, 1, 2, • • • . The sign ahead 11^ can be "+" or 
"-" when p e = 0. 

From the above discussion, we know that one root of equation M = is needed in 
the variable transformation, namely fi tpi . To avoid the complexity caused by introducing 
complex, we always use the real one. Luckily, equation B M = always has real roots, but 
which is not true for equation R(r) = 0. For cases in which equation R(r) = has no real 
roots we will use the Jacobi's elliptic functions sn(z|/c 2 ), cn(z\k 2 ) to express r. 




3.3. r coordinate 

If equation R(r) = has real roots, then r tpi exists, we can define the following constants 
by using r tpi : 

60 = 4r t 3 pi - 2(q + A 2 - a 2 )r tpi + 2[q + (A - a) 2 ], (35) 

b 1 = 2r 2 tpi -^(q + X 2 -a% (36) 
4 

b 2 = -r tpi , (37) 

h = 1, (38) 
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and introduce a new variable t, 



bo 1 61 

4(r-r tpi ) 4 



It is similar with /i, using t as the independent variable, we can reduce the r part of equation 
(24) into the standard form of Weierstrass' elliptical integral 

f t{r) dt 

P =± 77T, I » ( 4 °) 

J V 4t - 92t - 93 

where g 2 = \{b\ — b b 2 ), g% = j^(Sb bib 2 — 26f — 60&3). Taking the inverse of above equation, 
we get t = p{p ± IL-; g 2 , #3)- Solving equation (|39|) for r, we have 



= A 7 t rr 6 ° — \ — r + r *Pi ' ( 41 ) 

4p(p±n r ;c/ 2 ,c/3) - h 

where 

n r = |p -1 [*(r<m);^2,^3]|- (42) 

The sign ahead IT also depends on the initial value of p r , which is the r component of four 
momentum of a photon, and 

Pr > 0, +, 

Pr = 0, { r ™ = r ^ ^= ™ ' +'-' (43) 
[r ini = r tp2 , U r = (z+njuj , +, -, 

Pr < 0, 

where 00 is the period of p(z; g 2 , #3) and n — 0, 1, 2, • • • . The sign ahead II r can be "+" or 
" -" when p r = 0. 

If equation R(r) = has no real roots, we use the Jacobi's elliptic functions to express r. 
Since the coefficient of r 3 is zero, the roots of equation R(r) = satisfy r\ + r 2 + r^, + = 0. 
Therefore the roots 7*1, r 2 , r$, can be written as 

r\ = u — iw, r 2 = u + iw, 

r3 = — u — iv, Ti = —u + iv. (44) 
Introducing two constants Ai and A 2 

\ 1 2 = J—[4u 2 + v 2 + w 2 ± J{4u 2 + w 2 + v 2 ) 2 - 4w 2 v 2 }, (45) 
which satisfy Ai > 1 > X 2 > 0, and a new variable t, 



/ Ai - 1 / Ai + 1\ 

* = (Ax -AdKr + l r " "AT^lJ ' (46) 
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we can reduce the r part of equation (24) to the Legendre's standard elliptic integral 

ft 

P= ' 



= 1 



dt 



where 



Ai — A2 



m 2 



Ai 



(47) 



(48) 



—2u ± w(Xi 


— A 2 )sn(j9^V 


% ± n \m 2 ) \cn{pw^/\[ ± n \m 2 )\ 


(Ai 


— A 2 )sn 2 (]9K; 


V / A^±n |m 2 )-(A 1 -1) 



Using the definition of Jacobi's elliptic function sn(z|/c 2 ) ( Abramowitz fc Stegunfl965 ), from 
equation (47) we obtain t = sn(pw^/X^ ± n |m 2 ). Solving the equation (46) for r, we get the 
expression of r as the function of p 

r±(p) u 

where 



(49) 

n 2 ) - (Ai - i) 

n = [*(r i7li )|m 2 ] |. (50) 

When the initial value of p r > 0, we have r = r_, and when p r < 0, we have r = r+. And p r 
can not be zero, otherwise the initial radial coordinate of the photon will be one root of 
equation R(r) = 0, which is the case that has been discussed above. 



3.4. t and (f) coordinates and affine parameter a 

In this section, we will express the coordinates i, (f) and the affine parameter a as the 
numerical functions of the parameter p. All of these variables have been expressed as the 
integrals of r and 9 in the equations (|9])-( 11 ) and (17)-(22). The goal is achieved if we can 



compute all of these integrals along a geodesic for a specified p. Making transformations 



from r and fi to a new variable t (defined by equations (30), ( |39| ) and (|46|)), we will compute 
these integrals under the new variable t. For simplicity we use F r (t) and Fe(t) to denote the 
complicated integrands (see below) for r and 9 respectively. 

Firstly, we discuss the integral path, which starts from the initial position and terminates 
at the photon. If the photon encounters turning points along the geodesic, then the whole 
integral path is not monotonic, as shown in Figurejl] for poloidal motion (radial motion 
is similar). In this figure the projected poloidal motion of a photon onto the r-9 plane is 
illustrated. The motion is confined between two turning points: ji tpi and /Li tp2 . The photon 
encounters the turning points for three times. Obviously any sections of the path which 
contain one or more than one turning points is not monotonic, such as path CDE, EFP etc. 
The path between any two neighboring turning points has the maximum monotonic length 
and the total integrals should be computed along each of them and summed. 
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There are four important points involved in the limits of these integrals, i.e., ji tpi) \i tp2 , 
Hi n i and n P) they are the \i coordinates of turning points, initial point and the photon position 
for a given p respectively. And the values of these points corresponding to the new variable t 
are t tpi , t tp2) t ini , which can be calculated from equation (30), and t p) which can be calculated 
from t = p(p±n M ; g 2 , gs) with a given p. Because the function t = t{ji) expressed by equation 
(30) is monotonically increasing, we have t tp2 < t ini < t tpi and t tp2 <t p < t tpi . 

If we use Nti and Nt 2 to denote the number of times of a photon meeting the turning 
points ji tpi and fi tp2 for a given p respectively, and define the following integrals (cf. Figure 

I„= (" F,(t)dt, h = / F,(t)dt, I 2 = j F,(t)dt, 

hi = h + h= I F 9 (t)dt, I 02 = I 2 -I = F e (t)dt, (51) 



Up. 



•2 



then the integrals of 9 in a, t and (j) then can be written as (cf. Figure [T]) 

cj e (fa t e ) = -sign( Pe )I + 2Nt 1 I 1 + 2Nt 2 I 2l 

= -[sign(^) + 2Nh - 2Nt 2 ]I + 2Nt 1 I 01 + 2Nt 2 I 02 , (52) 

where p$ is the 9 component of the initial four momentum of a photon. In order to evaluate 
the above expression, we need to know Nt\ and Nt 2 for a given p. Similarly if we define five 
integrals from the equation (31) as follows: 

f tp dt f ttp i dt f tp dt 

Po = / /— Pi = / 1= -v , P2 



ini VW)' * J tp vW K 2 y/W®' 

pm ^ P0+p ^ L vwy ^^-"^ L TWr (53) 

where W(t) = 4:t 3 — g 2 t — g 3 , and we will get the following identity: 

p = -sign(p0)p o + 2iV£ipi + 2Nt 2 p 2) 

= -[sign(^) + 2Nh - 2Nt 2 ]p + 2Nt lPoi + 2Nt 2 p 02 . (54) 

And notice that Nti and Nt 2 are not arbitrary and related to the initial direction of the 
photon in poloidal motion. For p$ > (or p$ = and 0^ = 9 tpi ), they will increase as 



Nt x = 1 1 2 2 • • • , 
Nt 2 = 1 1 2 2 3 • • • . 
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For po < (or p$ = and ini = tp2 ), they will increase as 



Nt x = 1 1 2 2 3 • • • , 
iVt 2 = 1 1 2 2 • • • . 

For a given we find that there always exists one pair of Nt\ and Nt 2) which satisfy equation 
(54) and they are the number of a photon meeting the turning points. With Nt\ and Nt 2) 
the equation (52) now can be evaluated readily. 

For r coordinate, the process is similar with the above. Nt\ and Nt 2 also represent 
the number of times of a photon meeting the turning points r tpi and r tp2 respectively. Five 
integrals are defined as: 



Ioi 



F r (t)dt, h 

rtt Pl 



h + h 



F r (t)dt, I 2 = 
F r (t)dt, I 02 = I 2 -I 



F r (t)dt, 



F r (t)dt, 



(55) 



Up 2 



where t p = p(p ± II r ; g 2) g^) (or t p = sn(pwy/X^ ± n |m 2 ) when equation R(r) = has no 
real roots), and t piJ t P2 and t ini are calculated from equation (39) or (46). Then the integrals 
of r in a, i, cj) can be written as 



(3 ' y ( ijry 



) = -sign(p r )I + 2Nt 1 I 1 + 2Nt 2 I 2 , 

= -[sign(p P ) + 2^! - 2Nt 2 ]I + 2A^t 1 / i + 2M 2 1 02 . 



(56) 



To get Nti and Nt 2) we define p , pi and p 2 from equation (40) as: 



Po 



^ tint 



dt 



Pi 



/ 

Jtr 



Pi 



dt 



P2 



dt 



P01 =Po+Pi 
For a given p, we have 



tpi 



dt 



Vwttj' 



P02 =P2~P0 



ini dt 



ttp 2 



Vmtj' 



(57) 



p = -sigii(p r )p + 2Ntipi + 2Nt 2 p 2 , 

= - [sign(p r ) + 2NU - 2Nt 2 ]p + 2Nt lPoi + 2Nt 2 p 02 . 



(58) 



To get Nt\ and Nt 2 from above equation one just needs to notice that when p r > (or 
p r = and r ini = r tpi ) they will increase as 

Nt! = 1 1 2 2 • • • , 

7V£ 2 = 1 1 2 2 3 • • • , 
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when p r < (or p r = and r ini = r tp2 ) they will increase as 

Nh = 1 1 2 2 3 • • • , 
iV* 2 = 1 1 2 2 • • • . 



Similarly for a given p, there is one pair of Nt\ and Nt 2 satisfies equation (58). With Nt\ 



and Nt 2 equation (56) now can be evaluated. In many cases the number of a photon meeting 



the turning points in r is less than 2, especially when r tpi is less or equal to r^, or r tp2 is 
infinity, or equation R(r) = has no real roots, both Nt\ and Nt 2 will be zero. 



3.5. Reductions to Carlson's elliptic integrals 

In previous sections, four coordinates r, 0,0, t and the affine parameter a have been 
expressed as functions of p, and in which many elliptic integrals need to be calculated. In 
this section, we shall reduce these integrals into standard forms and then evaluate them by 
Carlson's method as Dexter & Agol ( 2009| ) did. 

Firstly, we introduce two notations Jk(h) and Ik(h) with following definitions: 

Jk{h) = / , (59) 

rr 2 f J r 

h{h) = / , (60) 

J ri (r — h) k yj (r 2 — 2ur + u 2 + w 2 )(r 2 + 2ur + u 2 + v 2 ) 

where k is an integer. From equation ([8]), we get one of the standard forms as 

Jo = / t= = • (61) 
J tl ^^t 6 - g 2 t - g 3 



After being reduced to Jo, the forms of integrals of r and 6 in (|8j) are exactly same. Noticing 
the definition of parameter p, we have J = P> The radial integrals in equation ^ are 
reduced to 

^ = ^ j2 uJ+^ Ji UJ +r ^ (62) 

where J has been replaced by p. The radial integrals in equation (10) can be reduced to 
U = <J r + | Ji (^) + (2r tpi + 4 + A t+ - A t _)p - B t+ .h(t + ) + Bt_J.it.), (63) 
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where 



r± = 1 ± y/l — o? ) 

. 2[r ± (4-aA)-2a 2 ] 
At± = 



B t ± = 



U 



(r+ - r-)(r tpi - r±) ' 

[r±(4-aA) -2a 2 ]6 
2(r+ - r_)(r tpi - r ± ) 2 ' 

61 ^0 



4 4(r±-r tpi ) 

Similarly the radial integrals in the equation ( [TT| ) have the following form 

(j) r = a[(A 0+ - A^p - B n Ji(t + ) + B^_Jx (<_)], 

where 

2r± — aA 



(64) 



(65) 



0± 



(r+ - r_)(r tpi - r±) ' 
(2r± — aA)6 



= 77 ~ ~ " ,<V " U r^. (66) 

4(r + - r_)(r tpi - r±) 2 

When equation R(r) = has no real roots, the integrals of r in a, t and 4> can be written as: 



where 



Or = 1-2(0), 

t r = a r + 4p + 2/_x(0) + C t+ h(r + ) - C t -h(r_), 
<f> r = a[C^ + /i(r+) - C^_/i(r_)], 



2[r±(4-aA) -2a 2 ] 2r± - aA 
w± — " " > ^4>± — • 



r + — r_ r + — r 

The integrals concerning /i in the equation ^ are reduced to 



a^ = a 



4 J + ^~ Jl 



7" 
^2 



16 



and ta = a a- The \i integrals in the equation (11) can be reduced to 



= A 



where 



-^ + W, + J 1 (t + )-W,_J 1 (t.) 

1 fl tp 



WW = 



(67) 
(68) 
(69) 

(70) 
(71) 
(72) 



t± = 4 + 



1 ±^pi) 2: 

&0 



4 4(±i-/i tpi y 



(73) 



- 15 - 



Finally we have 



t t r + t^, 



A 6 ' 



(74) 



From equations (61)-(72) we know that the integrals need to be calculated are J , Ji, J 2 , 
and Ii, I-i, 1-2- Now we use the Carlson's method to evaluate them. When equation 
4t 3 — g 2 t — 5^3=0 has three real roots denoted by e\, e 2 and e^, Jk(h) can be written as 



(Carlson 1988) 



i r 

Jk{h) = s h - j 



dt 



^/{t-e l ){t-e 2 ){t-e z ){t-hf k ' 



= ^-[-1,-1,-1,-2*], 



(75) 



where Sh=sign[(y — h) k ]. When equation 4t 3 — g 2 t — g$ = has one real root t\ and one pair 



of complex conjugate roots u + iv and u — iv, Jk(h) can be written as ( |Carlson 1991) 

dt 



Jk(h) 



Sh 



2/ 



2j y y/(t - e^it 2 - 2ut + u 2 + v 2 )(t - h) 



2k ■ 



s fc ^[-l,-l,-l,-2fc]. 



(76) 



From the equations (30) and (39), we know that when r = r tpi or \i = /x tpiJ t will be oo, thus 
one limit of these integrals can be infinity. 



The integrals Ik(h) can be reduced to Carlson's integrals directly ( |Carlson||1992[ ) 
Sh I 

J V 



h(h) 



dr 



y \J (r 2 — 2ur + u 2 + w 2 ){r 2 + 2ur + u 2 + v 2 ){r — h) 2k ' 
s h [-l,-l,-l,-l,-2fc], (77) 



where [pi, • • • ,Pk] is a symbol used by Carlson to denote the elliptic integrals and has the 
following definition: 



[Pw- M= f Y[(ai + bit)^ 2 dt, 

J v i=l 



(78) 



Carlson 


(1988, 


1989, 


1991, 


1992) 
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Constants of motion 



4.1. Basic equations 



In previous sections we have expressed all coordinates as functions of a parameter p and 
discussed how to calculate them by Carlson's method. But before the calculation one needs 
to specify the constants of motion and p r) P9) which determine the signs ahead n^II^IIo 
and also how the number of turning points increasing. In this section we shall discuss how 
to compute A and q and p r) pe from p( a ), which are the components of the four-momentum 
measured in the LNRF reference and have been specified by the user. 



Firstly, following Bardeen et al. (1972) we introduce the LNRF (locally nonrotating 
frame) observers or the ZAMO (zero angular momentum observer), the basis vectors of the 
orthonormal tetrad of them are given by 



e (a) (LNRF) = e^, 



(79) 



where 



-to 



V 








ue~ u 
















e"^ 2 








e"^ 



(80) 



/ 



And the covariant components of the four momentum of a photon in the B-L coordinate can 
be expressed as 



P„ = E[ -1, 



(81) 



where s r and sg are signs of r and 6 components. One can easily show that p( a ) 



namely (Shakura 1987) 



Pit) 

P(r) 

Pie) 
P(4>) 



-Ee 



A 



(1-Acu), 
R 



= s e Ee-» 2 - 
= EXe~^. 



(82) 

(83) 

(84) 
(85) 



From equations (83) and (84) we have 5 r =sign(j9( r )) and S0=sign(p(0)), which determine the 



initial direction of the photon in the B-L system (cf. equations (34) and (J43J) ) , thus determine 
the way how the number of the turning points increasing. 
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Solving equations (82) and (85) simultaneously for A, one obtains 

^ = sin ep(4>)/p(t) 



AE/A + io sin 9p w /p (t) 



(86) 



Using A and equation (82), one obtains E = p^e v /{l — Xco). Using A and E, from equation 

rula of i 



(84) one obtains the formula of calculating the motion constant q, 

2 



AE/A + usmdp^/p(t) 



a 



cos 2 



.Pit) 



(1 - Xco) 



A 



(87) 



Thus we have obtained the basic equations (86) and (87) connecting A, q and the components 
of four momentum p^ measured in the LNRF reference. When p^ are given the constants 
of motion and the initial direction of the photon are both uniquely determined. 

To prescribe p^ , one should notice that they satisfy following equation 



P(t) + P(r) + P(9) + P(<f>) 



0, 



thus there are only three independent components. Obviously the user can specify the 
four momentum p^ directly in LNRF or equivalently specify p',^ in anyother reference 
frame of his/her own choice and then to transform it to the LNRF reference by a Lorentz 



transformation, i.e., p( a ) = o^p'^y where cta J is the transformation matrix. From equations 
(86) and (87) we know that what one needs is just p(i)/p(t) an d 

«? ) + «?V w /Pit 



pa) 



\t) 



4 t} + ^>P' {J) /p{ t ) 

The cffl should be specified by the user according to his/her needs. ^ 

As an example, in Figure [2| we show a group of null geodesies emitted isotropically 
from a particle moving around a black hole in a marginally stable circular orbit (r ms ) with 



Si). 



(89) 



2 When a reference frame K' has physical velocities v t ,vq, with respect to a LNRF, the general Lorentz 
transformation matrix has six independent parameters, i.e., = a£\6i,02,03]v r ,VQ,v ( f ) ), where 6i are 
the angles between the corresponding spacial basis vectors of the two references. If 0i = 0, the matrix can 
be written as follows (Mis ner et al.||1973 ) 



7 ~l v r 

-~/V r l+7 2 ^ 2 /(l+ 7 ) 

-jv e 7 2 ^r/(l + 7) 

-7^0 7 2 ^r/(l + 7) 



-<yv e 
j 2 v r v e /(l + 7) 
l + 7% 2 /(l + 7 ) 

7 2 ^^/(! + 7) 



-7^ 
7 2 ^r*V(l + 7) 
7 2 ^^/(! + 7) 
l + 7 2 ^/(l + 7) J 



(90) 
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a=0.9375. And the physical velocities of the particle with respect to the LNRF are v T — vq — 



— 



of the partic 



-uj 



The four-momentum p'^ are specified isotropically in the reference 



e and then transformed to the LNRF by the Lorentz transformation expressed 
by equation (90), i.e., ^( a ) = da P{b)- Withp( a ) the constants of motion are computed readily. 



The light bending and beaming effects are illustrated obviously in this figure. 

In the next section, we shall discuss how to compute P(i)/p(t) from impact parameters, 
which play a key role in imaging. And for simplicity we shall use the transformation expressed 



by the equations (90) and (91) if the observer has motion 



4.2. Calculation of motion constants from impact parameters 



From the works of Cunningham & Bardeen (1973) and Cunningham (1975), we know 



that A and q can be calculated from impact parameters, usually denoted by a, /?, which are 
the coordinates of the hitting position of a photon on the photographic plate of the observer. 



The formulae provided by them read as follows (Cunningham & Bardeen 1973) 



A -asin0 ofe5 , 

q = f3 2 + (a 2 — a 2 ) cos 2 O & 5 . 



(92) 
(93) 



The above equations are valid only when the distance between the observer and the emitter 
is infinite and the observer is stationary. Practically the distance is not infinite, otherwise the 
integrals of coordinate t will be divergent. When the distance is finite, the above formulae 
should be modified. We extend those formulae to general situations, in which both the finite 
distance and the motion state of the observer are considered. 

To consider the finite distance is very easy. One just needs to substitute the coordinates 



Tobs^obs of the observer into equations (86) and (87) in the calculation. While to consider 



the motion state is more complicated. Obviously we can distinguish the motion states of the 
observer into two kinds. In the first one the observer is stationary and in the second one the 
observer has physical velocities v r ^v$^v$ with respect to the LNRF reference. 



where 7 = [1 — (1% 



A)] x / 2 , and its inverse form 



a 



(a) 



( 7 7^r 

JVq 7 2 ^r/(l + 7) 

\ 7^0 7 2 ^r/(l +7) 



JV0 

<y 2 v r v e /(l + 7) 
l + 7 %7(l + 7) 
7 2 ^*V(! + 7) 



7^0 \ 
7VV(l + 7 ) 

7 2 wW(l + 7) 
l + 7 2 ^/(l + 7) J 



(91) 
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In the first kind, the observer is just a LNRF observer and whose orthonormal tetrad is 
given by equation (79), namely e( a )(obs) = e( a )(LNRF). While in the second kind, the tetrad 
of the observer can be created by a Lorentz transformation, i.e., e( a )(obs) = S^e^LNRF). 
Here 5^ is given by the equation (91). 

As shown in Figure |3| we plot the image of a photon hitting on the photographic plate. 
The plate is located in the plane determined by the basis vectors e^)(obs) and e(^)(obs), and 
in which an orthonormal coordinate system a, /3 has been established. The basis vectors e a 
and ep of the system are aligned with e(^)(obs) and e^)(obs) respectively. All photons will 
go through the center of the Lens before hitting on the plate. From this figure, one obtains 
the relationships between the impact parameters and p',^ as follows: 



a rscal 



P'(r) 



(3 rscal 



P{r) 



r=r obs ,t 



r=r obs , 



(94) 
(95) 



And obviously one can read off p'^ > 0, ap'^ > and fip'^ > 0. 



Two dimensionless factors r and seal have been multiplied to amplify the size of the 
image, otherwise which will be infinite small, since the distance D between the central 
compact object and the observer and the size of the target object L satisfy D 3> L. 



In the rest frame of the observer, the spacetime is locally flat, we still have 

-P(t)+P(r)+P(0)+PW =0- 



(96) 



Using equations (94), (95) and (96) and noting the signs of p 7 ^, a, /3, we obtain 



P'(r) _ 


1 


P{t) 


Jl + (a/r seal) 2 + (/3/r seal) 2 


P{o) _ 


/3 /rscal 


P'(t) 


\Jl + (a/r seal) 2 + (/3/r seal) 2 


p'w _ 


a /rscal 


p'(t) 


y 1 + (a/r seal) 2 + (/3/r seal) 2 



(97) 
(98) 
(99) 



Substituting equations (97)-(99) into (89), we get the functions /3)/p(t)(a, f3). We can 

then calculate A and q from impact parameters by using equations (l86j) and (87). 
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One can verify directly that when r 5 s oo, seal = 1, and v r = vq — v§ — 0, the 



equations (86) and (87) reduce to (92) and (93) immediately. 



If the observer has motion, the image on the plate will have a displacement compare 
to the image when the observer is stationary. The displacement is proportional to the 
observer's velocity and can be described by a C) /3 C , which is the coordinates of image point of 
the origin of B-L coordinate system on the photographic plate. Obviously a c and /3 C satisfy 
the following equations: 



P(0)(a c ,p c ) = 0, 
P(0)(a C3 /3 c ) = 0. 



(100) 
(101) 



Actually, P(</>)(a, /3) = represents the projection of the spin axis of the black hole onto the 
plate. When v r — vq — v^ = 0, the above equations become a c = and f3 c = 0. The region 
of the image on the plate therefore is [—AL + a C) AL + a c ] and [— AL + /3 C , AL + /3 C ], where 
AL is the half length of the image. 



4.3. Redshift formula 



The redshift g of a photon is defined by g = E b s /E ern . From above discussion, we know 
that E b s = —p[ t y E em = —p^u^ m) where u^ m is the four- velocity of the emitter. If we define 



It — OL t + a t — . 

P(t) 



(102) 



From equation (82), we have — p' t = Ee u (l — Xoo)/ f t . Using the equation (81), then g can 
be expressed as follows 



9 



[e-"(l - Xu)/f t ] 



obs 



(103) 



u* (l + s r rVR/A + seBy/We - Afi) 

where f = u r /u l ,0 = u° /u t ) Vt = (j) = u^/u 1 are the coordinate velocities. 

With r,#, fi, the physical velocities of the emitter with respect to the LNRF v r ,VQ,v$ 



can be written as (Bardeen et al. 1972) 



r, v Q = e^e, v^ = e^~ u (n - u) 



with which the four-velocity of the emitter can be expressed as 



(104) 



(105) 
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Then the g can be rewritten as ( Muller fc Camenzind||2004 ) 

[e-"(l - Au;)// t ] o6s 



5 



ye-" (1 + s r e v v r \/R/y/TA + s e e u v e ^/W e /\^ - XQ 

For an emitter movies in a Keplerian orbit, the formula of g reduces to 

[e-"(l - Au;)// t ] o6s 



(106) 



5 



[ 7 e-" (1 - XQ)l 



(107) 



5. A brief introduction to the code 



5.1. The four coordinates and affine parameter functions 

We have expressed the four coordinates r, /i, 0, t and the affine parameters a as functions 
of p. We denote them as follows: 

r(p),Kp), </>(p),t(p),a(p). (108) 

In practical applications, we are interested in determining the original position where the 
photon was emitted or the regions traveled by the photon. To make the calculations effec- 
tively, all photons are traced backward from the observer to the emitter along the geodesies. 
But not all photons start from the observer will go through the emission region one inter- 
ested, and the tracing process will be terminated either these photons go to the infinity or 
fall into the event horizon of a black hole. 

Now we discuss how to determine the intersection of a geodesic with the surface of an 
optically thick emission region, we assuming that the optical depth of which is so large that 
a sharply emission surface exits. And the surface is smooth and continuous and can be 
described by an algebra equation: 

F(r,0,0) = F o , or F{x,y,z)=F , (109) 

where x ) y, z are the pseudo Cartesian coordinates and defined by 

x = \Jr 2 + a 2 sin 6 cos 0, y = Vr 2 + a 2 sin 6 sin 0, z = rcos9. (110) 

In some special cases the surface one considered may not keep stationary, the surface equation 
will be a function of time £, i.e., F(r, 6 ) 0, t) = F . We introduce a function f(p) defined by 



f(p) = F[r(p),0(p), ( f>(p)]-F o . 



(Ill) 
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Then the roots of equation f(p) — correspond to the intersections of the geodesic with 
the target surface. Therefore if equation f(p) = has no roots, the geodesic will never 
intersect with the surface. To solve this equation effectively, we classify geodesies into four 
classes denoted by A, B, C and D, according to their relationships with respect to a shell, 
shown in Figure [5j The shell includes the emission region completely, and its inner and outer 
radius are r in and r out . Reminding that r tpi and r tp2 are turning points, between which the 
radial motion of a photon is confined. Geodesies in the four classes satisfy the conditions A: 
r tpi > r ouU B: r in < r tpi < r ouU C: r h < r tpi < r in and D: r tpi < r h respectively, where r h is 
the radius of the event horizon. 

The values of parameter p corresponding to the intersections of the geodesic with the 
shell are denoted by pi ) p 2 , p^ ) p^ which satisfy p\ < p 2 < P3 < Pa- Obviously the roots of 
equation f(p)=0 may exist on intervals [pi,P2] an d [P3,P4\- We use the Bisection or the 
Newton- Raphson method to search the roots ( Press et al.||2007 ). 



Solving the radiative transfer equation in optically thin or thick media, one needs to 
evaluate integrations along geodesic with taking the affine parameter a as the independent 
variable. Since we have taken p to be the independent variable, we can replace a by p to 



evaluate these integrals. From the definition of i.e., equation (23), one has 



dr d6 

p= 7W) = <112) 

and from the equations ([I]) and ^ one gets 

da = ±E^L= = ±E^L (113) 

From above equations one immediately obtains 

da = Srfp, (114) 
which converts the independent variable from a to p in radiative transfer applications ( |Yuan 



et al. 2009). 



Finally we give a brief discussion on the determination of a geodesic connecting the 
emitter and observer ( Viergutz||1993 ; Beckwith fc Done]|2005 ) . We use a emi ft em to represent 



the impact parameters of the geodesic which connecting the observer and emitter and p em to 
indicate the position of the emitter on the geodesic, in which the coordinates of the emitter 
and (j) em . Obviously we have the following set of equations: 

(116) 

(f)(p 

em) ^emi Pem) 0em- (H7) 
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In principle if we can solve this set of nonlinear equations simultaneously for p emi a em) /3 em , 
the geodesic is determined uniquely. Therefore the observer-emitter problem also becomes a 
root finding problem. In our code we use the Newton- Raphson method (Press et al. 2007) 
to solve these equations. 



5.2. The code 



In this section we shall give a brief introduction for the code, and a more detailed intro- 
duction is given in the README^] file. The code is named YNOGK (Yun-Nan Observatory 
Geodesies Kerr) and written by Fortran 95, in which the object-oriented method has been 
used. The code is composed by a couple of modules. For each module, a special function has 
been implemented and one can use all supporting functions and subroutines in that mod- 
ule by a command "use module-name" in his/her own program. By adding corresponding 
modules into one's own code, one can easily develop new ones to handle special and more 
sophisticated applications. 

Two modules named ell-function and BLcoordinate are the most important ones 
in ynogk, the former one includes supporting functions and subroutines for calculating 
the Carlson's elliptic integrals and the R-functions. Many routines in this module come 



from geokerr (|Dexter fc Agol 2009) and Numerical recipes (Press et al. 2007). The latter 



module includes routines for computing all coordinates and the affine parameter functions: 
t{p) ) r{p) ) 6{p) ) (j){p) and cr(p). To call these routines, constants of motion and components 
of four-momentum p( r ), p^, p^ measured in a LNRF reference must be prescribed, which 
can be computed by two subroutines named lambdaq and initialdirection in ynogk. The 
former routine computes the constants of motion from impact parameters, while the latter 
one computes A and q from the initial p f ^ given in a reference K\ which has physical veloc- 
ities Vr^vo^Vcfy with respect to the LNRF. Of course one can compute them by his/her own 



subroutines according to their needs. According to the discussion in Section [5T| we present 
a module named pern- finding to search the minimum root p em of equation f(p) = 0, and 
f(p) as an external function should be given by the user. In module obs-emitter, we present 



routines to find the root of equations ( 115 )-( 117) by using the Newton-Raphson algorithm 



( Press et al.|2007 ). In the testing section of the code, this module has been used to determine 
the geodesic connecting the observer and the central point of a hot spot, which movies in 
the inner most stable circular orbit (ISCO). With the geodesic the motion of the spot can 
be described easily. The results are agree very well with previous works, in which a very 



^http : / / wwwl . ynao . ac . cn/~ yangxl/ readme . pdf 
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different method has been used to determine the motion of the spot, i.e., by tabulating the 



motion according to the time of the observer over one period (Schnittman & Bertschinger 



2004; Dexter fc Agol 2009). 



All routines for computing the Carlson's elliptic integrals have been extensively checked 
by NIntegrate function of Mathematica. The original code of these routines comes from 



geokerr (Dexter & Agol 2009), and has been modified to adapt to our code. The original 



code for the computing of R-functions comes from Numerical recipes (Press et al. 2007). 



The same check also has been done for the functions t(p), r(p), 0(p), 0(p), o"(p). When \a\ 
and |/3 1 < 10 -7 , we let them to be zero, since the Carlson's integrals can not maintain their 
accuracy. The treatment is same for any other parameters if they take offending values. For 
some critical cases, special treatments also have been implemented. 



5.3. Comparisons and speed tests 



Our code has many common points with geokerr of Dexter & Agol (2009). We both 



use the Carlson's method to compute the elliptic integrals and use the elliptic functions to 
express all coordinates as functions of a independent variable. But the elliptic functions we 
used are mainly the Weierstrass' elliptic function p{z\ #2, #3), which has a cubic polynomial, 
leading simpler root distribution and the cases of integral are reduced. In our code the four 
B-L coordinates r, 0, 0, t and affine parameter a are expressed as analytical or numerical 



functions of a parameter p : which corresponds to I u or in Dexter & Agol (2009). With 



this treatment one can compute the geodesies directly without providing any information 
about the turning points in advance. Which also allows one to track emissions from a more 
sophisticated surface, not only for standard thin accretion disk. In the code testing section 
we will show the images of a warped disk, which has a curved surface. 

Our strategy, i.e., expressing coordinates as functions of p semi-analytically, can be 
extended to compute the timelike geodesies directly, almost without any modifications. As 



mentioned in Dexter & Agol (2009), the calculations of the timelike geodesies involve many 



more cases. The main challenge is to specify the number of radial turning points for bounded 
orbits in advance. But our strategy does not require the specification of the number of turning 
points both in radial and poloidal coordinates in advance, therefore which can be used 
naturally and effectively in the calculations of timelike geodesies even in a Kerr-Newmann 
spacetime. 



In our code we give the orthenormal tetrad of the emitter or the observer analytically, 
provided the physical velocities of which with respect to the LNRF are specified. They 
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may be useful in Monte-Carlo type code of radiative transfer, which needs one to make 
transformations from the reference of the emitter to the B-L coordinate system frequently 
(Dolence et al. 2009). As illustrated in figure [5J emissions in the reference of the emitter 



are specified isotropically, but from the perspective of the B-L coordinate system which are 
anisotropic due to the Doppler beaming effect. 



Our testing results for various toy problems agree well with those of Dexter fc Agol 



(2009). In Figure |4| we illustrate the projection of a uniform orthonormal grid from the 
photographic plate of the observer onto the equatorial plane of a black hole, in which the 
solid and dotted lines represent the results from our code and geokerr respectively. They 
agree with each other very well. 

The basic strategy used in ynogk to compute the elliptic integrals and functions are 
very similar to geokerr. For example we make ynogk to compute the minimum number of R- 
functions possible and share them between routines. This strategy improves the speed of our 
code greatly. But there is still some differences between the two codes. Firstly, we assemble 
r(p) and ji{p) into routines for computing t{p) ) (j)(p) and cr(p), thus the repeated calculations 
for the same integrals among those functions can be avoided. We provide a routine named 
ynogk to compute the four B-L coordinates and affine parameter simultaneously. We also 
provide two independent routines named radius and mucos to compute r(p) and n(p) 
respectively. Secondly, ynogk can save the values of variables used in the calculations for a 
same geodesic but for different p. These values can be used repeatedly. 

ynogk has almost same speed with the geokerr in tracing radiations from an optically 
and geometrical thin disk, because there is only one point, i.e., the intersection of the ray 
with the disk surface, needs to be calculated. For the calculations of radiative transfer, in 
which many points are needed along each geodesies, ynogk has a little slower than geokerr. 
The speed tests of a code are not only dependent on the applications mostly, but also on 
the environment of the code running. From the testing results of ynogk, we expect that the 
speed of which is almost same with geokerr in many other applications. 

For more detailed introductions, one can see the README file. In the next section, we 
will show the testing results of our code for toy problems. 



6. The tests of our code 

In this section, we shall present the testing results of our code for toy problems. The 
results not only demonstrate the validation our code, but also give specific examples of its 
utility. Firstly, we show the image of a black hole shadow, in which the intensity represents 
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the value of the affine parameter a. By this example we want to test the validation of 
function a(p). Then we show the images of a couple of accretion disks and a rotationally 
supported torus, all of them are optically thick and have a sharply emission surface. The 
disks include the standard thin, thick and warped disks. Next we show the images of a ball 
orbiting around a Kerr black hole in a Keplerian orbit to illustrate the gravitational lensing 
effect. Then we calculate the line profiles of the Fe Ka and the blackbody radiation spectra 
of a standard thin accretion disk around a Kerr black hole. In order to test the accuracy of 
function t(p), we image a hot spot moving around a Kerr black hole in the ISCO for various 
black hole spins. We also calculate the spectrogram and light curves of the spot over one 
period of the motion with various inclinations for a Schwarzschild black hole. Finally we 
discuss the radiative transfer equation and its solution, whit which the radiative transfer 
process in a radiation dominated torus around a black hole has been discussed. We give 
the images of the torus for optically thin and thick cases. The resolution of images in this 
section is taken to be 801 x 801, each pixel corresponds to an unique geodesic. 



6.1. Black hole shadow 

As the first test of our code we give the image of a black hole shadow. We trace 
all photons backward from the photographic plate to the black hole along geodesies. The 
intensities of the image are taken to be the affine parameter a evaluated from the observer to 
terminations on the geodesic — either when it intersects with the event horizon of the black 
hole or reaches a turning point and returns to the starting radius. The evaluations of the 
affine parameter outside the shadow are multiplied a factor 1/2. In Figure [6j we show the 
image from an edge-on view. We take the spin a to be 0.998, and the distance of the observer 
to be 10 6 r g , where r g is the gravitational radius. 

To evaluate affine parameter a from function o"(p), we need which is the value of 
parameter p corresponds to the event horizon and also the root of equation r(p) = rv We 
can get ph by evaluating the integral of r in the definition of p ) and need not to solve this 
equation directly. We provide a routine named r2p to complete this evaluation. 



6.2. Accretion disks 



Next we present the images of the accretion disks around a Kerr black hole, including 
the standard thin, thick and warped disks. The imaging of the disks is usually taken as 
the first step to calculate the line profiles of the Fe Ka and the spectrum (Li et al. 2005). 
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Usually the pseudo colors of the image represent the redshift g or the observed flux intensity 
l v of emissions come from the disk. 

In Figure [7j we show the image of a standard thin disk, the inner and outer radius of 
which is r ms and 22 r g repectively. The black hole spin a is 0.998 and the inclination angle 
Oobs is 86°. The distance of the observer is 40 r g . The shape of the image is quite different 
from the one observed from infinite far away. We also illustrate the high-order images of the 
disk in this figure. Due to the light bending and focusing, one can see the part behind the 
black hole and the bottom side of the disk. The color intensities represent the redshift g of 
emissions come from the disk. 

In order determine the intersections of geodesies with the disk, we need to know the min- 
imum root of equation fi(p) = 0. We provide two routines named pemdisk and pemdisk-all 
to compute the root by evaluating the integral of \± in the definition of p. Using pemdisk one 
can draw the direct image, while using pemdisk-all one can draw the direct and high-order 
images. 

The surface of the thick disk has a constant inclination angle S with respect to the 



equatorial plane ( |Wu fc Wang 2007). To trace the thick disk, we need to solve equation 
ji(p) = cos(7r/2 — 5) to get p for the upper surface and ji(p) = cos(7r/2 + 5) for the bottom 
surface, the roots of these two equations can also be computed by pemdisk and pemdisk- 
all. Since the surface particles of the disk no longer keep in the equatorial plane, they will 



do the sub-Keplerian motion with a angular velocity given by (Ruszkowski & Fabian 2000) 



tt/2 



i - 



tt/2 



(118) 



where VLk = l/(r 3 / 2 + a) is the Keplerian velocity and parameter n is taken to be 3 here. 
Using the equation ( 107[ ), we can calculate the redshift of the emissions come from the disk. 



The images are shown in Figure pi which agrees very well with Figure 10 of Wu fc Wang 



(2007) 



The warped accretion disk is also a very interesting object in astrophysics (Bardeen & 



Petterson 


1972; 


Wu & Wang 


2007; 


Wang & Li 


2012) 



for the warped disk, in which the disk is assumed to be optically thick and its surface can 



be described by (Wang k Li 2012) 



cot 9 = — tan f3 cos(0 — 7) , 



(119) 
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where parameters 7 and /3 are defined by 

l{p) 70 + niexp n 2 - 
7r r(p) 



— r(p) 



n 3 sin 



^ out T'in 
Tin 



(120) 



(121) 



. 2 T out fir, 

where r in and r out are the inner and outer radius of the disk, and ni, n 2 , n 3 are the warping 
parameters. With above equations, we get the f(p) as follows 



f(p) = tan P(p) cos[(f)(p) - 7(p)] 



(122) 



>/l-0 2 (p)' 

With the minimum root of equation f(p) = 0, we can image the warped disk. For the 
poloidal velocity 9 of the particle is nonzero, the formula (103) or (106) is used to calculate 
the redshift g (cf. [Wang fc Li| ( [20121 ). The images of the warped disk are shown in Figure 
[9j in which the warping parameters n\ and n 2 are nonzero, leading the disk warps along 
azimuthal direction. For comparison one can see Figure 3 of Wang & Li (2012), in which 774 
is taken to be zero for simplicity, thus the shape of the disk is quite different from the one 
illustrated here. 



6.3. Rotationally supported torus 



In this section, we give the images of a rotationally supported torus. For simplicity we 
give a brief introduction for the torus model here, for the more detailed discussions one is 



recommended to the paper of Fuerst & Wu (2004) or Younsi et al. (2012). The torus is 



assumed to be stationary and axisymmetric. Due to the balance of the centrifugal force, 
gravity and pressure force, the structure of the torus is stratified and the isobaric surfaces 



where 



can be described by a set of differential equations (Younsi et al. 2012) 

dr ip2 

d6 _ -^1 



= M 



E - 2r 2 



E 2 
Mr, 



^ 2 = sin 26 [ ^[afT 1 - (r 2 + a 2 )] 2 + 



asin#) +rsin 2 #, 
A" 



ft = 



E 2 
M 



[r sin 6>) 3/2 + a\[M V r sin 6 



(123) 
(124) 

(125) 
(126) 
(127) 
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and Q is the angular velocity. r& represents the radius at which the the particle orbits with 
a Keplerian velocity. The index parameter n is crucial for regulating the angular velocity 
profile and adjusting the geometrical aspect ratio of the torus. ( is an auxiliary parameter. 
In order to give the outer surface of the torus, one needs to specify the most inner radius 
of the torus, which is usually regarded as the intersection of the isobaric surface with the 
ISCO. Taking the inner most radius in the equatorial plane to be the initial condition, the 
differential equations (123) and (124) are now readily to be integrated. In Figure 10 we 



illustrate the images of the torus, which has the same parameters with Figure 3 of Younsi 



et al. (2012), the results agree with each other very well. 



6.4. The gravitational lensing effect 

Due to the strong gravity field, when the trajectory of a photon is closed to the vicinity 
of a compact object, it will be bent or focused, then multiple images will be observed, this 
is the so called gravitational lensing effect. Here this phenomenon will be illustrated by 
a simple example, in which a ball moves around a near extremal black hole (a=0.998) in 
a Keplerian orbit. The radius of the orbit is i? , then the angular velocity of the ball is 
Cl = l/(i?^ 2 + tt). The coordinates of the center of the ball will be 



x Q (p) = \ Rl + a 2 cos[ttt(p)], (128) 



y (p) = ^R 2 + a^m[nt(p)i (129) 
*b(p) = 0. (130) 

Then the function of the surface of the ball can be expressed as follows: 

f(p) = Vi x (p) - %o(p)} 2 + Hp) - vo(p)} 2 + Av? - (131) 

where r\ is the radius of the ball. The images observed from an edge-on view are illustrated 



in Figure 11 For different positions of the ball in its orbit, the image changes greatly, which 



even becomes a ring as the ball movies to the back of the even horizon. 



6.5. The line profiles of Fe Kcr 

The calculation of line profiles is very easy provided the structure of the disk is specified. 
For simplicity, we assume that the particles of the accretion flow do the Keplerian motion 
and the disk is geometrical thin and optically thick. The inner and outer radius of the disk 
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are located at r ms and 15 r g respectively. The emission is monochromatic and the profile of 
which can be described by the Dirac's 8 function in the local rest frame of the flow 



5{v 



(132) 



where n is the index of emissivity and assumed to be 3. Since l v jv 3 is an invariance along a 
geodesic ( Misner et al.||l973 ), we get the observed intensity / 5 S = < g ,3 / em , where g = v bsl v era 
is the redshift. Then the observed flux density F v at frequency v can be computed by 
integrating 7 5 S over the whole plate as following expression 



e 



$(v ~ v e m)g 3 dadf3. 



(133) 



The observed intensities have been normalized in the computation. The results are shown 
in Figure 12, which agrees very well with the Figure 3 of Cadez et al. (1998). From this 



figure one can see that the higher black hole spin leads the broadening in low frequency for 
the ISCO is closer to the event horizon, the gravitational redshift effect is remarkable, while 
the higher inclination angle leads the broadening in high frequency for the Doppler beaming 
effect. 



6.6. The blackbody radiation spectrum of a Keplerian disk 



In this section we will compute the spectrum of a Keplerian disk around a Kerr black 
hole to illustrate effects of the black hole spin and the observer's inclination angles on the 



observed profiles of the spectrum (Li et al. 2005). Similarly, the disk is assumed to be 



geometrical thin and optically thick, and the radiation spectrum of the disk in its local rest 
frame is an isotropic blackbody spectrum. We denote the effective temperature of the disk 
by T e ff. Then radiation intensity at frequency v can be written as 



Iemiy) 



his 3 



exp(hv / k B T e $) - 1' 



(134) 



where h and k B are the Plank and Boltzmann constants respectively. For a blackbody 
radiation, the effective temperature is simply 



^eff 



F(r) 

°SB 



1/4 



(135) 



where <jsb is the Stefan-Boltzmann constant. Here we do not consider effect of the returning 
radiation of the disk on the spectrum, therefore F(r) is just the energy flux emitted from 
the disk's surface measured by a locally corotating observer. For the Keplerian accretion 
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disk around a Kerr black hole, Page & Thorne (1974) have get the analytical expression for 
F{r), i.e., 



M 



(136) 



where M is mass accretion rate, / is a function of r, a,r ms , and the seminal expression 



of which is given by equations (15d) and (15n) of Page & Thorne (1974). With F(r) the 
effective temperature of the disk can be computed readily. Using the invariance /^/^ 3 , one 
can get the observed intensity 7 5 S = g 3 iem- The total observed flux density at frequency v 
therefore is the integration of l Q \) S over the whole plate 



F{v bs) 



I 



exp(hv obs /gk B T eS ) - 1 



dad/3. 



(137) 



Then the photon number flux density is AT & 5 = F(yob 8 )/ Eob 8 . 



The results are plotted in Figure 13 Compare to the Figure 5 of Li et al. (2005) we 



find that the basic features of the two figures are in agreement. For example, as shown in 
top panel, we see that as the spin of the black hole goes up the spectrum becomes harder. 
Physically, this is due to the fact that as the spin a increases, the system of the accretion disk 
has a higher radiation efficiency and a higher temperature. In the bottom panel of Figure 
131 we can see that at the low-energy end, the flux density goes down as # 5 s goes up. As 



explained by Li et al. (2005) this is caused by the projection effect. While at the high-energy 



end, the flux density goes up as the O 5 5 increases. As pointed out by Li et al. (2005) this is 



resulted from the joint action of the effects of Doppler beaming and gravitational focusing. 

In the top panel there is a noticeable effect: even though we do not consider the returning 
radiation, the flux density goes up as the spin increases in the low-energy end. |Li et al. 



(2005) suggested that this effect is caused by the returning radiation. We proposal that this 



explanation may be not correct and the effect is just caused by the simple fact that a higher 
spin leads to a higher radiation efficiency and temperature. 



6.7. The motion of a hot spot 

In order to test the validation of the function t(p) and illustrate the time delay effect 
in the Kerr spacetime, we image a hot spot orbiting around a black hole retrogradely in a 
marginally stable circular orbit for various spins and compute the observed light curve and 
spectra. The radius of the spot is i? sp ot=0.5 r g . The emissivity of the spot is taken to be 
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Gaussian shape in its rest frame (Schnittman & Bertschinger 2004), i.e., 



j(x) oc exp 



x — X, 



■spot 



(*)l s 



2R 



spot 



(138) 



For the motion of the spot, one must consider the time delay effect and the azimuthal position 
of the spot when imaging the spot and calculating its spectra. In order to compute the time 
delay At for each geodesic starting from the photographic plate, a reference time i b s , which 
is taken to be the time used by a photon traveling from the central point of the spot to the 
observer, needs to be specified. Meanwhile the position of the spot can be determined by 
its central coordinates (r ms , /x = 0, sp ot)- Then with the method discussed in section [5j 
(i.e., the method to determine a geodesic connecting the observer and emitter with the given 
coordinates) , we can determine the geodesic connecting the central point of the spot and the 
observer. With this geodesic the reference time £ b s can be calculated readily. Using £ b s , 
we can easily calculate the time delay At for each geodesic, and At = t geo — i b s , where £ geo 
is the time used by a photon traveling from the observer to the disk following the geodesic. 
With At and the position of the spot, we can compute the distance between the intersection 
of the geodesic with the disk and the center of the spot, i.e., |x — x spot |. Thus the emissivity 
can be computed readily. 



In Figure [HJ we illustrate the images of the spot with different black hole spins. As 
the spin increases, the marginally stable circular orbit is closer to the event horizon of the 
black hole, and the time delay effect becomes remarkable. The image of the spot is seriously 
warped, especially when the spot movies to the back of the event horizon. 

When an image is obtained, the redshift and Gaussian emissivity of all points on the 
spot can be computed. Repeating this procedure over one period of the motion gives a 
time-dependent spectrum. Integrating the spectrum over frequency, or equivalently over the 
impact parameters, gives the light curve. The spectrum and light curve are shown in Figure 



15, which agree well with the results shown in Figures 6 and 7 of Dexter & Agol (2009). 



6.8. Radiative transfer 



6.8.1. The radiative transfer formulation 
In this section we give a brief discussion to the radiative transfer process under the Kerr 



spacetime. One can find more detailed discussions from Fuerst & Wu (2004) and Younsi 
et al. (2012). It is well known that X = I^/is 3 , x = VOL v an d V = i v jv 2 are Lorentz invariants, 



where l v is the specific intensity of the radiation, a v and j v are the absorption and emission 
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coefficients at the frequency v. The radiative transfer equation reads (Younsi et al. 2012) 

dl 



— = -1+^. 



(139) 



where r v is the optical depth at the frequency v, and defined by dr v = a v ds and ds = 
—p^u^da, in which ds is the differential distance element of a photon traveling in the rest 
frame of the medium, a is the affine parameter, p^ is the four momentum of the photon, and 
■u M is the four velocity of the medium. Then the radiative transfer equation can be rewritten 
as (lYounsi et al.||2012|) 



(140) 



The solution of above equation is ( | Younsi et al. 2012) 



1(a) = l(a )e- T ^ - 
where the optical depth is 



[ 3 ~^P- ex P (" [„ MOIpXM^') P^U„da", (141) 



a. v (a')p tl u' 1 \ (J ida' . 
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^0 



As discussed in section [5TT| we can convert the independent variable from affine parameter a 
to parameter p. Using a = a(p) and da = Srfp, we can rewrite the solution as the integration 
of parameter p (Yuan et al.| 2009) 



X(p) = l(p )e- T ^ 
where 



p Mp") 

^^exp 



PO 



a v (p')\pX\v'\Z'dv' pX\p"^"dp", (143) 



a v (p')PpU ,1 \ p iE'dp'. 



(144) 



PO 



With above formulae one can deal with radiative transfer problems without considering the 
scattering contributions to the absorption and emission coefficients as did by [Yuan et al. 



(2009 


) and 


Younsi et al. 


(2012) 



6.8.2. Radiative transfer in pressure supported torus 



In section |6.3| we have discussed a rotationally supported torus and demonstrated its 
images. When the torus is optically thick, only the emissions come from the boundary surface 
are considered. When the torus is optically thin, all parts of the torus will do contributions 
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to the observed emissions. We need to consider the radiative transfer procedure along the 
ray inside the torus. To get the absorption and emission coefficients, we need to konw the 
structure model of the tours, which determines the distributions of the temperature, mass 
density, pressure etc. 

Firstly we construct the model of the torus, in which the torus is a perfect fluid and its 
energy-momentum tensor is given by QYounsi et al. 2012) 



T a ? = (p + P + e)u a u p + Pg a ^ 



(145) 



where p is the mass density, P is the pressure, and e is the internal energy, u a is the four 
velocity of the fluid, and g a/3 are the contravariant components of the Kerr metric. From the 
conservation law, namely T a ^ = 0, we get the equation of motion of the fluid as follows 
(lAbramowicz et al.||1978|): 



d n P 



-U a .j3U h 



(146) 



p + P + e 

where the semicolon ; represents the covariant derivative, and u a tfvP 
acceleration of the fluid. For the torus is stationary and axisymmetric, we have a t = 0, 



a n is the four 



0, and a r , a# are given by (lYounsi et al. 2012) 
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where i = it* is the time component of the four- velocity, f2 is the angular velocity and takes 
the form of equation (127). They satisfy following equation 

1 



u 



y/-{9tt + + 



(149) 



Since the torus is assumed to be radiation dominated, the pressure P can be regarded as the 
sum of gas pressure P gas and radiation pressure P ra d, and 



pk B T 



gas 



P ra.d 



pmu 
aT 4 



0P, 



= (1 - P)P, 



(150) 
(151) 



where &b is the Boltzmann constant, \i is the mean molecular weight, mn is the mass of a 
hydrogen, (3 is the ratio of gas pressure to the total pressure, and a = 7T 2 k 4 /15h 3 c 3 is the 



black-body emission constant. From the above equations, one finally obtains 

45(1-/5) "" 1/3 



p = hc 



kT = he 



7r 2 (^m H /5) 4 
45(1-/5)"" 1/3 



P 



4/3 



7T 2 /XmH/5 



p 



1/3 



(152) 
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Thus P = np T \ which implies that the state equation of the fluid is polytropic, therefore its 
internal energy is proportional to the pressure e — P/(T — 1), and the equation of motion of 



the fluid (146) becomes 



P + 



P ) a n = —d ry P. 



r- 1 

Substituting d a P = /^rp r_1 9 a p and P = Kp v into above equation, one obtains 

^2-r 



(154) 



d a p = -a a 



P 



P 



(155) 



kT r - i . 

Introducing a new variable ^ defined by £ = ln(r — 1 + nYp 1 * -1 ), above equation is simplified 



as 



(156) 



which implies that the vector n = (a r) ao) in the r-9 plane can be regarded as the normal 
vector of the contours of the density p. Thus if we use t = (rfr, d6) to denote the tangent 
vector of the contours, we have n • t = 0, or equivalent ly 

a r dr + aedO = 0. (157) 

If we use ds to denote the differential proper length of the tangent vector, we have 

ds 2 = t • t = g rr dr 2 + g ee d6 2 , (158) 

where g rr and gee are the components of the Kerr metric, and g rr — E/A, gee — S. Solving 
the equations (157) and (158) simultaneously, we get a set of differential equations to describe 



the contours of density p 
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If we introduce an auxiliary variable ( defined by d( = ^ l /~KfEds, and substitute equations 



(147) and (148) into the above equations we get 



dr 
~d~C 
d6 
d( 



(161) 
(162) 



which have the exactly same forms with equations (123) and (124), where tpi and ip2 are 



given by the equations (125) and (126). With these equations, the distributions of the mass 



density p of the torus now are readily to be computed by evaluating the integral of £ from 
the torus center r = r^, p = p c to the location (r, 6) along a path C which is orthogonal to 
the density contours everywhere. And the integral of £ is 



a r dr + aedO. 



(163) 



c 



From the equations (152) and (153) one can get the total pressure and temperature distri- 



butions immediately with the given density p. 

Knowing the structure model of the torus, the absorption and emission coefficients are 
now readily to be specified, with which we can discuss the radiative transfer process inside 
the torus. Using the above torus model, we shall give two examples of radiative transfer 
applications. 

Firstly we consider a rather simple case, in which the torus is optically thin. The 
emissivity is taken to be proportional to the mass density p, i.e., j em oc p, and is independent 
on the frequency v. The absorption coefficient a v is simply assumed to be zero. The 
torus parameters are n = 0.21, = 12 r r The black hole spin a is 0.998. The ratio of 



gas pressure to total pressure /3 is 2.87 x 10 -8 . In Figure 16 we draw the images of the 
torus, which is optically thin and radiation pressure dominated. We see that the emissions 
mainly come from the central region of the torus, where the density is higher. As the 
inclination angle of observer increases, the frequency shift of the emission caused by the 
Doppler boosting becomes larger. In this figure the false color represents the observed 
intensities of the emission, showing that the approaching side of the torus is brighter than 
the receding side especially at higher inclination angles. 

Secondly, we mimic a more realistic case, namely the thermal free- free emission and 
absorbtion procedure, in which the torus is semi-opacity. The emission and absorbtion 
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coefncients of the torus for a photon at energy E are given by (Younsi et al. 2012) 



a(Eo) 



= K 



Bi 



cm~' 
cm~ 
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-1/2 



keV J 



(164) 
(165) 



where B = k&T, /C and B\ are the normalization constants, n e is the electron number 
density and n e = p//xmn, (Tt is the Thompson cross-section. The observed intensity images 



of the optically thick and semi-opacity torus are plotted in Figure [17} These images are 
quite different from those of an optically thin torus. The emissivity now depends on the 
temperature, which decreases towards to the outer surface of the torus, leading the limb 
darkening phenomenon. When the rays are nearly tangential to the layers of the torus, they 
will travel a longer distance and go through the outer, thus colder layers. While when the 
rays are perpendicular to the layers of the torus, they will travel a shorter distance and 
go through the inner therefore hotter layers. Consequently, the observed intensity at lower 



inclination angles will be much brighter than that at higher inclination angles (Younsi et al. 
2012b. 



7. Discussions and conclusions 

Following Dexter & Agol (2009|) we have presented a new public code named ynogk for 



the fast calculating of null geodesies in a Kerr spacetime. The code is written by Fortran 95, 
and composed by a couple of modules. In which the object-oriented method has been used, 
which makes the addition of the code to one's own readily. 

In ynogk the B-L coordinates r and \± have been expressed as analytical functions of the 
parameter p. In these expressions, the Weierstrass' and Jacobi's elliptic function p[z\ #2, #3)5 
sn(z|/c 2 ) and cn(z|k 2 ) are used, since the reductions to Weierstrass's standard integrals are 
much easier, in which only one real root of the equations R(r) = and B M = is required. 
The B-L coordinates £, (j) and the affine parameter a have been expressed as numerical 
functions of p. For a given p, the number of times of a photon reaches the turning points 
both in radial and poloidal motions is uniquely determined and needs not to be specified by 
the user. 

Actually in addition to p, coordinates r, /1 (or 9) can also be taken as the independent 
variables ( [Dexter fc Agol|2009l ) . The main reason of using p is that one can pay no attention 



to handle turning points, which has been done by the inner routines of our code. This virtue 
is convenient for a person who is not familiar with or has no interesting to the details of 
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the calculation of a geodesies in the Kerr spacetime. Another reason is that the value of 
p which corresponds to the termination of the geodesic — either at the infinity or the event 
horizon — is finite. Thus it is easier to handle p than r. In our code r and \i can also be taken 
as the independent variable. We provide a routine named geokerr, which can take r or \i 
as the independent variable. But the number of turning points should be prescribed. 

With the expressions of all coordinates and affine parameter as functions of the ray- 
tracing problem, which determines the intersection of the ray with a target object, now 
becomes a root finding problem. The function f(p) that describes the surface of the target 
object needs to be given by the user and the roots of equation f(p) = correspond to the 
intersections. We provide a module named pem-finding to search the minimum root of this 
equation by the Bisection or the Newton-Raphson method. In addition, the observer-emitter 
problem can also be converted to a root finding problem, which requires one to solve a set of 
nonlinear equations. A module named obs-emitter based on the Newton-Raphson method 
to solve these equations is provided in our code. The routines in this module will return the 
solution, provided the coordinates of the emitter, r em) 6 ern and em , are given. 

We present a new set of formulae to compute the constants of motion A and q from 
initial conditions. These formulae can be regarded as the extensions of |Cunningham fc 



Bardeen (1973). Our formulae are pervasive and can be used to handle more sophisticated 
cases, in which the motion state and the finite distance of the observer or the emitter with 
respect to the black hole are considered. One may find it is convenient when dealing with 
problems in which the emitter has motion and is closed to the vicinity of a black hole, e.g., 
the self-irradiation process in the inner region of a disk. 

The code has been tested extensively with various toy problems in the literature. The 
results agree well with previous works. The comparisons with geokerr of |Dexter fc AgoT 



(2009) also have been presented. 



Finally we point out that the strategy discussed in this paper can be naturally extended 
to the calculation of the timelike geodesies almost without any modification. Especially for 
the timelike bounded orbits, in which the number of turning points both in poloidal and 
radial coordinates can be arbitrary. The extension of this strategy to calculate the timelike 
geodesies in a Kerr-Newmann spacetime has been done and the results are under preparation. 
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Fig. 1. — This figure illustrates the motion of a photon in the 9 coordinate, which has been 
projected onto the r — 9 plane. The motion is confined between two turning points fi p i 
and fi P 2. A, D and F indicate the positions where the photon reaches the turning points 
and P indicates the position of the photon. The path between any two neighboring turning 
points (such as DA, DF) has the maximum monotonic length and the integrals of 9 should 
be evaluated along each monotonic section and summed. The doted (such as CA, EF and 
PF) and solid (such as DC and DE) lines represent the integral paths of Ii and I 2 (see text) 
respectively. Obviously the BC section is the integral path of and one has I = 



J i = Ic = Ie = Jp > J 2 = Id = Id> etc - 
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Fig. 2. — A set of geodesies emitted isotropically from a particle orbits around a black hole in 
the marginally stable circular orbit with a=0.9375. x and y are pseudo-Cartesian coordinates 
in the equatorial plane of the black hole. The figure shows the light bending and beaming 
effects clearly. A circle in the center represents the boundary of the event horizon. 
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Fig. 3. — The figure shows the hitting of a photon on the photographic plate of the observer, 
from which the relationships between the impact parameters a, (3 and the components p'^ 
of the four momentum of the photon are derived. Before hitting the plate, all photons will 
go through the center of the lens. e( r )(obs), e^)(obs) and e(^)(obs) are the contravariant 
basis vectors of the frame, and the basis vectors of a, /3 coordinates e a , are aligned with 
e^)(obs), e(0)(obs) respectively. 
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Fig. 4. — The projection of a uniform grid from the photographic plate of the observer onto 
the equatorial plane of a black hole is shown. The inclination angle 9 ohs is 60° and the black 
hole spin a is 0.95. Solid lines represent the results from our code and the dotted lines from 
geokerr. x and y are pseudo-Cartesian coordinates in the equatorial plane of the black hole. 
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Fig. 5. — The classifications of a set of null geodesies according to their relationships with 
respect to a shell, the inner and outer radius of which are r in and r out respectively. The 
geodesies are classified into four classes, marked by A, B, C and D. Since the target object 
or the emission region are assumed to be completely included by the shell, only geodesies 
in classes B, C and D have probabilities to intersect with the target object or go through 
emission region. The trajectories of the geodesies are schematically plotted and have been 
projected onto the r-9 plane. The central black region represents the black hole shadow. 
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Fig. 6. — The shadow of a black hole with near extremal spin (a = 0.998) from the edge-on 
view is shown. The radial coordinate of the observer is 10 6 r g . The greyscale represents the 
value of the affine parameter a evaluated from the observer to the terminated position — 



either at the black hole or re-emerging to the starting radius. Compare to Figure 2 of Dexter 



& Agol (2009). a and /3 are the impact parameters, which describe the size and the position 



of the image on the photographic plate. 
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Fig. 7. — This figure shows the image of a standard thin accretion disk, whose inner and 
outer radius are r ms and 22 r g respectively. The black hole spin a is 0.998 and the inclination 
angle O 5 S is 86°. The radial coordinate of the observer is 40 r g . One can see Figure 6 of 



Beckwith & Done (2005) or Figure 4 of Dexter & Agol (2009) for comparison. The high-order 



image is also shown, a and f3 are the impact parameters, and the intensity of the greyscale 
represents the redshift g of emissions come from the surface of the disk. 



- 49 - 




Fig. 8. — Images of a thick disk around a near extremal Kerr black hole (a=0.998) for 
various inclination angles are shown. The surface of the disk has a constant inclination angle 
5 with respect to the equatorial plane and 5 is taken to be 30°. The inner radius is the 
marginally stable circular orbit r ms and the outer radius is 20 r g . The inclination angles 
Oobs are 5°, 30°, 55° and 80° for panels a, b, c and d respectively, a and ft are the impact 
parameters, and the intensities of the color represent the redshift g of emissions come from 
the surface of the disk. Compare to Figure 10 of |Wu fc Wang (2007). 
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Fig. 9. — This figure shows the images of a warped accretion disk around a near extremal 
black hole (a=0.998) viewed from different azimuthal angles. The inner and outer radius 
of the disk are r ms and 50 r g . The observer's inclination angle # 5 S is 50°. The warping 
parameters are rt\ = 47r, n 2 = 4, and n 3 = 0.95. The azimuthal angle 70, which represents 
the view angle, is 0°, 45°, 90°, 135°, 180°, 225°, 270° and 315° for panels from left to right 
and top to bottom respectively. For comparison we show a image observed from a face-on 
view in the final panel. The false color also represents the redshift g of the emissions come 
from the surface of the disk, a and /3 are the impact parameters. We take the parameter 
rii 7^ 0, leading the warping of the disk along the azimuthal direction shown clearly in the 



final panel, which is the main difference compare to Figure 3 of Wang & Li (2012). 
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Fig. 10. — Images of a rotationally supported torus, which is geometrical and optically 
thick, are shown. The torus parameters are n = 0.2, r& = 12 r g . The black hole spin a is 0, 
0.5 and 0.998 for panels from top to bottom. The inclination angle O 5 S is 45° for left column 
and 85° for right column. The false color represents the redshift g of the emissions come 
from the surface of the torus and the white areas represent the zero-shift regions, a and /3 
are the impact parameters. Compare to Figure 3 of Younsi et al. (2012). 
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Fig. 11. — This figure illustrates the gravitational lensing effect by the motion of a ball 
movies around a near extremal black hole (a=0.998) in a Keplerian orbit. The motion 
is observed from an edge-on view. The radius of the ball and the orbit are 5 r g and 20 
r g respectively. The central red region represents the black hole shadow. The azimuthal 
angles of the ball measured from the line of sight along inverse-clockwise direction are 
0°, 90°, 150°, 160°, 180°, 195°, 210°, 270° for panels a-h. The pseudo color shows the redshift 
of the emissions come from the surface of the ball, a and /3 are the impact parameters. 
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Fig. 12. — The theoretical line profiles of the Fe Ka of a thin accretion disk for various black 
hole spins and inclination angles are shown. The inner and outer radius of the disk are r ms 
and 15 r g . Top row: a = 0.2; bottom row; a = 0.998. Left column: # 5 S = 10°; middle 
column: O 5 5 = 30°; right column: # 5 s = 75°. The horizonal and vertical axes represent 
the frequency and flux of the line respectively and are normalized. The line profiles are in 
agreement well with the Figure 3 of Cadez et al. (1998). 
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Fig. 13. — The effects of the black hole spin (top) and the inclination angle (bottom) on the 
spectrum of a standard thin accretion disk around a Kerr black hole are shown. The inner 
and outer radius of the disk are r ms and 30 r g , where r ms is the radius of marginally stable 
circular orbit. 
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Fig. 14. — The images of a hot spot orbits around a black hole for different black hole spins 
are shown. The spot lies in a standard thin accretion disk and its central point is fixed at 
the ISCO. We have extended the inner radius of the disk to the photon orbit r^, at which 
the energy per unit rest mass of a particle is infinite. It is also the innermost boundary of 
circular orbit for particles ( Bardeen et al.fl972 |). For the panels from left to right and top to 
bottom, the black hole spin a is 0.998, 0.5, and -0.998 respectively. The inclination angle 
Oobs is 85°. The false color represents the value of g — j(x), where g is the redshift of the 
emissions come from the surface of the disk, and j(x) is the emissivity of the spot. 




Fig. 15. — The time-dependent spectrogram (panel a) and light curves (panel b) of a hot 
spot orbits around a Schwarzschild black hole in the marginally stable circular orbit (6 r g ) 
over one period are shown. The inclination angle o b s is 60° for the spectrum. The greyscale 
in panel a represents total sum of emissivity j(x) of emissions which are observed at the same 
time and have the same redshift g. The greyscale has been normalized and the maximum is 
taken to be 1. Compare to Figure 6 and 7 of Dexter & Agol (2009). 
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Fig. 16. — This figure shows the images of an optically thin and radiation pressure dominated 
torus. The inclination angles of the observer are 15°, 30°, 45°, 60°, 75° and 90° for panels 
from left to right and top to bottom. The black hole spin a is 0.998, and the ratio of gas 
pressure to total pressure f3 is 2.87 x 10 -8 . The torus parameters are n = 0.21, = 12 r g . 
The brightness of each pixel represents the observed intensity integrated along a geodesic 
ray at a given frequency and has been normalized, and the maximum for each panel is the 
same and equals to 100. a and /3 are the impact parameters. 




Fig. 17. — This figure shows the images of an optically thick and semi-opacity torus. The 
inclination angles are 15°, 30°, 45°, 60°, 75° and 90° for panels from left to right and top 
to bottom. The black hole spin a is 0.998, and the ratio of gas pressure to total pressure 
/3 is 2.87 x 10 -8 . The torus parameters are n = 0.21, r& = 12 r g . The brightness of each 
pixel represents the observed intensity integrated over the entire spectrum. The intensity 
has been normalized, and the maximum of each panel is the same and equals to 100. a and 
/3 are the impact parameters. 



